require(phyloseq)
require(tidyverse)
require(phyloseq)
require(reshape2)
Loading required package: reshape2
Attaching package: ‘reshape2’
The following objects are masked from ‘package:data.table’:
dcast, melt
The following object is masked from ‘package:tidyr’:
smiths
require(dplyr)
require(ggplot2)
Load data
ps_dmn <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/16S/DMN_ests_16S.Rdata")
sample_data(ps_dmn)$Herbicide <- factor(sample_data(ps_dmn)$Herbicide, levels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
ps_rare <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/16S/HerbPt1_rare_16S.Rdata")
sample_data(ps_rare)$Herbicide <- factor(sample_data(ps_rare)$Herbicide, levels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
ps_trans <- readRDS("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/data/PhyloseqObjects/16S/HerbPt1_hel_trans_16S.Rdata")
sample_data(ps_trans)$Herbicide <- factor(sample_data(ps_trans)$Herbicide, levels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
create alphadiversity tables
alpha_div <- estimate_richness(physeq = ps_rare, measures = c("Observed", "Shannon", "Chao1"))
#pull out metadata and concatonate with alpha diversity metrics
md<-data.frame(sample_data(ps_rare))
alpha_div_md <- rownames_to_column(alpha_div, "Barcode_ID_G") %>% full_join(md)
Joining, by = "Barcode_ID_G"
alpha_div_md$Herbicide <- factor(alpha_div_md$Herbicide, levels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
Shannon Div plots - no significant differences among herbicide treatments at any of the three time points
ggplot(data = alpha_div_md, aes(Herbicide, Shannon, color= Herbicide)) + facet_grid(. ~ Time) + geom_boxplot() + theme_classic() + theme(axis.text.x = element_text(angle = 45, hjust = 1) )
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_Shannon.pdf")
Saving 7.29 x 4.51 in image
aov_t1<-aov(Shannon ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T1",])
plot(aov_t1$residuals)
summary(aov_t1)
Df Sum Sq Mean Sq F value Pr(>F)
Herbicide 4 0.0689 0.01723 1.345 0.266
Residuals 51 0.6533 0.01281
aov_t2<-aov(Shannon ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T2",])
plot(aov_t2$residuals)
summary(aov_t2)
Df Sum Sq Mean Sq F value Pr(>F)
Herbicide 4 0.0449 0.01123 0.865 0.491
Residuals 49 0.6355 0.01297
aov_t3<-aov(Shannon ~ Herbicide, data = alpha_div_md[alpha_div_md$Time == "T3",])
plot(aov_t3$residuals)
summary(aov_t3)
Df Sum Sq Mean Sq F value Pr(>F)
Herbicide 4 0.0457 0.01142 0.681 0.609
Residuals 50 0.8391 0.01678
ordinations and adonis testing with three separate objects (i.e., dmn, rarefied, transformed). Rare taxa are removed from rarefied and transfomred to sucessfully ordinate. At this point, the transformed data will not ordinate.
ord_dmn<-ordinate(physeq = ps_dmn, method = "NMDS", distance = "bray", k=3, trymax= 100)
ps_rare_sub<-prune_taxa(taxa_sums(ps_rare) > 2, ps_rare)
ord_rare<-ordinate(physeq = ps_rare_sub, method = "NMDS", distance = "bray", k=3, trymax= 100, previous.best = ord_rare)
#can't get the hellinger transformed data to ordinate successfully
ps_trans_sub<-prune_taxa(taxa_sums(ps_trans) > 0.5, ps_trans)
ord_transformed<-ordinate(physeq = ps_trans_sub, method = "NMDS", distance = "bray", k=3, trymax= 100)
Adonis testing of herbicide treatments by time point
ps_adonis<-function(physeq){
otu_tab<-data.frame(phyloseq::otu_table(physeq))
md_tab<-data.frame(phyloseq::sample_data(physeq))
if(taxa_are_rows(physeq)== T){
physeq_dist<-parallelDist::parDist(as.matrix(t(otu_tab)), method = "bray")}
else{physeq_dist<-parallelDist::parDist(as.matrix(otu_tab), method = "bray")}
print(anova(vegan::betadisper(physeq_dist, md_tab$Herbicide)))
vegan::adonis(physeq_dist ~ Herbicide * Time, data = md_tab)
}
ps_adonis(ps_rare_sub)
Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 4 0.002281 0.00057033 0.4454 0.7756
Residuals 160 0.204866 0.00128041
Call:
vegan::adonis(formula = physeq_dist ~ Herbicide * Time, data = md_tab)
Permutation: free
Number of permutations: 999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
Herbicide 4 0.7552 0.18880 1.1017 0.02657 0.037 *
Time 2 0.3792 0.18959 1.1063 0.01334 0.060 .
Herbicide:Time 8 1.5797 0.19747 1.1523 0.05559 0.003 **
Residuals 150 25.7059 0.17137 0.90450
Total 164 28.4200 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
ps_adonis(ps_trans_sub)
Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 4 0.01133 0.0028322 0.3169 0.8664
Residuals 165 1.47460 0.0089369
Call:
vegan::adonis(formula = physeq_dist ~ Herbicide * Time, data = md_tab)
Permutation: free
Number of permutations: 999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
Herbicide 4 0.2416 0.060388 0.97118 0.02277 0.573
Time 2 0.1771 0.088564 1.42433 0.01670 0.012 *
Herbicide:Time 8 0.5526 0.069078 1.11094 0.05209 0.067 .
Residuals 155 9.6378 0.062179 0.90845
Total 169 10.6091 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
ps_adonis(ps_dmn)
Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 4 0.019549 0.0048873 2.8407 0.02596 *
Residuals 165 0.283879 0.0017205
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Call:
vegan::adonis(formula = physeq_dist ~ Herbicide * Time, data = md_tab)
Permutation: free
Number of permutations: 999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
Herbicide 4 0.2226 0.055638 2.2900 0.04823 0.001 ***
Time 2 0.1044 0.052199 2.1484 0.02263 0.008 **
Herbicide:Time 8 0.5214 0.065169 2.6823 0.11299 0.001 ***
Residuals 155 3.7659 0.024296 0.81615
Total 169 4.6142 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Ordination plots DMN
ord_t1_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T1"), method = "NMDS", distance = "bray", k=3, trymax= 100)
Run 0 stress 0.0860376
Run 1 stress 0.08603687
... New best solution
... Procrustes: rmse 0.0006267826 max resid 0.003492555
... Similar to previous best
Run 2 stress 0.08603727
... Procrustes: rmse 0.000519486 max resid 0.002879424
... Similar to previous best
Run 3 stress 0.08603775
... Procrustes: rmse 0.0006710548 max resid 0.00365014
... Similar to previous best
Run 4 stress 0.08674028
Run 5 stress 0.08674051
Run 6 stress 0.08610309
... Procrustes: rmse 0.0155372 max resid 0.1010944
Run 7 stress 0.08662927
Run 8 stress 0.08674034
Run 9 stress 0.08603713
... Procrustes: rmse 0.0004570747 max resid 0.002500815
... Similar to previous best
Run 10 stress 0.1023472
Run 11 stress 0.08674014
Run 12 stress 0.09941847
Run 13 stress 0.08610274
... Procrustes: rmse 0.01547852 max resid 0.1009111
Run 14 stress 0.1017673
Run 15 stress 0.08662931
Run 16 stress 0.08610275
... Procrustes: rmse 0.01547835 max resid 0.1009027
Run 17 stress 0.08662977
Run 18 stress 0.08662964
Run 19 stress 0.08603738
... Procrustes: rmse 0.0005091999 max resid 0.002807653
... Similar to previous best
Run 20 stress 0.08610325
... Procrustes: rmse 0.01555667 max resid 0.1011618
*** Solution reached
T1_dmn<-ggordiplots::gg_ordiplot(ord = ord_t1_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T1_dmn$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_dmn_T1.pdf")
Saving 7.29 x 4.51 in image
ord_t2_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T2"), method = "NMDS", distance = "bray", k=3, trymax= 100)
Run 0 stress 0.1104159
Run 1 stress 0.1104149
... New best solution
... Procrustes: rmse 0.001171513 max resid 0.007161514
... Similar to previous best
Run 2 stress 0.1104149
... Procrustes: rmse 0.0002875201 max resid 0.001118723
... Similar to previous best
Run 3 stress 0.109526
... New best solution
... Procrustes: rmse 0.03334389 max resid 0.2369979
Run 4 stress 0.1095359
... Procrustes: rmse 0.0009442239 max resid 0.004563832
... Similar to previous best
Run 5 stress 0.1121088
Run 6 stress 0.1095264
... Procrustes: rmse 0.0001397644 max resid 0.0007317152
... Similar to previous best
Run 7 stress 0.1122111
Run 8 stress 0.1095263
... Procrustes: rmse 0.0008604593 max resid 0.003402465
... Similar to previous best
Run 9 stress 0.1095262
... Procrustes: rmse 0.0008348408 max resid 0.003298233
... Similar to previous best
Run 10 stress 0.1095262
... Procrustes: rmse 0.0005129767 max resid 0.00314258
... Similar to previous best
Run 11 stress 0.1108775
Run 12 stress 0.1247908
Run 13 stress 0.1104148
Run 14 stress 0.1095266
... Procrustes: rmse 0.0009426628 max resid 0.004072489
... Similar to previous best
Run 15 stress 0.1122114
Run 16 stress 0.1122109
Run 17 stress 0.1104146
Run 18 stress 0.1120888
Run 19 stress 0.1120856
Run 20 stress 0.1108786
*** Solution reached
T2_dmn<-ggordiplots::gg_ordiplot(ord = ord_t2_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T2")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T2_dmn$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_dmn_T2.pdf")
Saving 7.29 x 4.51 in image
ord_t3_dmn<-ordinate(physeq = subset_samples(ps_dmn, Time=="T3"), method = "NMDS", distance = "bray", k=3, trymax= 100)
Run 0 stress 0.09696447
Run 1 stress 0.09982855
Run 2 stress 0.09784799
Run 3 stress 0.09800983
Run 4 stress 0.09784894
Run 5 stress 0.0969634
... New best solution
... Procrustes: rmse 0.04061371 max resid 0.2635214
Run 6 stress 0.09706969
... Procrustes: rmse 0.04032592 max resid 0.2678095
Run 7 stress 0.09696485
... Procrustes: rmse 0.04056564 max resid 0.2625302
Run 8 stress 0.0969324
... New best solution
... Procrustes: rmse 0.004791775 max resid 0.03220976
Run 9 stress 0.09696451
... Procrustes: rmse 0.03977923 max resid 0.263354
Run 10 stress 0.09696596
... Procrustes: rmse 0.004560318 max resid 0.02412528
Run 11 stress 0.09784882
Run 12 stress 0.09838599
Run 13 stress 0.1005549
Run 14 stress 0.09777523
Run 15 stress 0.0969219
... New best solution
... Procrustes: rmse 0.001739056 max resid 0.008015476
... Similar to previous best
Run 16 stress 0.1008367
Run 17 stress 0.09697109
... Procrustes: rmse 0.03948731 max resid 0.2623541
Run 18 stress 0.09696141
... Procrustes: rmse 0.005146037 max resid 0.03412768
Run 19 stress 0.1028377
Run 20 stress 0.09707004
... Procrustes: rmse 0.03882214 max resid 0.2583721
*** Solution reached
T3_dmn<-ggordiplots::gg_ordiplot(ord = ord_t3_dmn, groups = data.frame(sample_data(subset_samples(ps_dmn, Time == "T3")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T3_dmn$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_dmn_T3.pdf")
Saving 7.29 x 4.51 in image
Ordination plots rarefied
ord_t1_rare<-ordinate(physeq = subset_samples(ps_rare, Time=="T1"), method = "NMDS", distance = "bray", k=3, trymax= 100)
Square root transformation
Wisconsin double standardization
Run 0 stress 0.1810491
Run 1 stress 0.181593
Run 2 stress 0.191872
Run 3 stress 0.1815709
Run 4 stress 0.1820694
Run 5 stress 0.1815727
Run 6 stress 0.1814585
... Procrustes: rmse 0.0482874 max resid 0.155698
Run 7 stress 0.1831936
Run 8 stress 0.1821565
Run 9 stress 0.1811472
... Procrustes: rmse 0.02328168 max resid 0.1241755
Run 10 stress 0.1809422
... New best solution
... Procrustes: rmse 0.01257791 max resid 0.07384735
Run 11 stress 0.1830541
Run 12 stress 0.1821389
Run 13 stress 0.1819163
Run 14 stress 0.1813678
... Procrustes: rmse 0.01983823 max resid 0.1171532
Run 15 stress 0.1821569
Run 16 stress 0.1830532
Run 17 stress 0.1819142
Run 18 stress 0.1837675
Run 19 stress 0.1860978
Run 20 stress 0.1811476
... Procrustes: rmse 0.01835014 max resid 0.1232744
Run 21 stress 0.1826836
Run 22 stress 0.182507
Run 23 stress 0.1817998
Run 24 stress 0.1825066
Run 25 stress 0.1815422
Run 26 stress 0.1819123
Run 27 stress 0.1816038
Run 28 stress 0.181148
... Procrustes: rmse 0.0183826 max resid 0.1233536
Run 29 stress 0.1915283
Run 30 stress 0.1823893
Run 31 stress 0.1821972
Run 32 stress 0.1827478
Run 33 stress 0.1837639
Run 34 stress 0.1823302
Run 35 stress 0.1856634
Run 36 stress 0.1837639
Run 37 stress 0.1811488
... Procrustes: rmse 0.01895138 max resid 0.1248699
Run 38 stress 0.180967
... Procrustes: rmse 0.02557771 max resid 0.1286128
Run 39 stress 0.183194
Run 40 stress 0.1837639
Run 41 stress 0.181913
Run 42 stress 0.1825652
Run 43 stress 0.1842313
Run 44 stress 0.1825074
Run 45 stress 0.1819123
Run 46 stress 0.1830533
Run 47 stress 0.1810718
... Procrustes: rmse 0.02863214 max resid 0.1313275
Run 48 stress 0.1917654
Run 49 stress 0.1831936
Run 50 stress 0.1834936
Run 51 stress 0.192229
Run 52 stress 0.182276
Run 53 stress 0.1822762
Run 54 stress 0.1825078
Run 55 stress 0.1815184
Run 56 stress 0.1821958
Run 57 stress 0.1814159
... Procrustes: rmse 0.04113413 max resid 0.1452096
Run 58 stress 0.1814153
... Procrustes: rmse 0.04101125 max resid 0.145104
Run 59 stress 0.186166
Run 60 stress 0.1814141
... Procrustes: rmse 0.0407537 max resid 0.1447673
Run 61 stress 0.181331
... Procrustes: rmse 0.03588519 max resid 0.1397534
Run 62 stress 0.1822761
Run 63 stress 0.1808178
... New best solution
... Procrustes: rmse 0.01546488 max resid 0.09409532
Run 64 stress 0.1811484
... Procrustes: rmse 0.02374035 max resid 0.1223041
Run 65 stress 0.1822021
Run 66 stress 0.1822765
Run 67 stress 0.1814159
Run 68 stress 0.1813262
Run 69 stress 0.1825039
Run 70 stress 0.1816037
Run 71 stress 0.1816037
Run 72 stress 0.1823672
Run 73 stress 0.182367
Run 74 stress 0.1810519
... Procrustes: rmse 0.02345899 max resid 0.1072637
Run 75 stress 0.181525
Run 76 stress 0.1825929
Run 77 stress 0.1810721
... Procrustes: rmse 0.03648172 max resid 0.1453981
Run 78 stress 0.1813256
Run 79 stress 0.1880286
Run 80 stress 0.1874007
Run 81 stress 0.185664
Run 82 stress 0.1814029
Run 83 stress 0.1808177
... New best solution
... Procrustes: rmse 0.000342235 max resid 0.001823038
... Similar to previous best
*** Solution reached
T1_rare<-ggordiplots::gg_ordiplot(ord = ord_t1_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T1")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T1_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_rare_T1.pdf")
Saving 7.29 x 4.51 in image
ord_t2_rare<-ordinate(physeq = subset_samples(ps_rare, Time=="T2"), method = "NMDS", distance = "bray", k=3, trymax= 100)
Square root transformation
Wisconsin double standardization
Run 0 stress 0.1534021
Run 1 stress 0.1512024
... New best solution
... Procrustes: rmse 0.0754476 max resid 0.2676264
Run 2 stress 0.1493358
... New best solution
... Procrustes: rmse 0.06201399 max resid 0.2214984
Run 3 stress 0.1532387
Run 4 stress 0.1493338
... New best solution
... Procrustes: rmse 0.03039888 max resid 0.1722904
Run 5 stress 0.1505376
Run 6 stress 0.1508442
Run 7 stress 0.1538361
Run 8 stress 0.1493389
... Procrustes: rmse 0.02035211 max resid 0.1193448
Run 9 stress 0.1493331
... New best solution
... Procrustes: rmse 0.03006402 max resid 0.1666307
Run 10 stress 0.1527413
Run 11 stress 0.1494084
... Procrustes: rmse 0.02890853 max resid 0.1471386
Run 12 stress 0.1503514
Run 13 stress 0.1564968
Run 14 stress 0.1493168
... New best solution
... Procrustes: rmse 0.01404071 max resid 0.07782122
Run 15 stress 0.1493021
... New best solution
... Procrustes: rmse 0.01741756 max resid 0.09981203
Run 16 stress 0.1493198
... Procrustes: rmse 0.02215659 max resid 0.1213156
Run 17 stress 0.1507341
Run 18 stress 0.149334
... Procrustes: rmse 0.00663877 max resid 0.03659659
Run 19 stress 0.1498004
... Procrustes: rmse 0.03768561 max resid 0.1683317
Run 20 stress 0.1495337
... Procrustes: rmse 0.03069104 max resid 0.1677622
Run 21 stress 0.1505776
Run 22 stress 0.1493306
... Procrustes: rmse 0.02690526 max resid 0.1437716
Run 23 stress 0.1516113
Run 24 stress 0.1498424
Run 25 stress 0.1493193
... Procrustes: rmse 0.02155654 max resid 0.1168637
Run 26 stress 0.1508459
Run 27 stress 0.1497011
... Procrustes: rmse 0.03123585 max resid 0.1580844
Run 28 stress 0.1499895
Run 29 stress 0.1494448
... Procrustes: rmse 0.007164553 max resid 0.03016139
Run 30 stress 0.1498771
Run 31 stress 0.1502641
Run 32 stress 0.1558884
Run 33 stress 0.1501174
Run 34 stress 0.1498656
Run 35 stress 0.1508896
Run 36 stress 0.1520205
Run 37 stress 0.1495492
... Procrustes: rmse 0.01320733 max resid 0.04627731
Run 38 stress 0.1493166
... Procrustes: rmse 0.002633167 max resid 0.01254371
Run 39 stress 0.1499219
Run 40 stress 0.1554524
Run 41 stress 0.1508795
Run 42 stress 0.1529847
Run 43 stress 0.1501535
Run 44 stress 0.1496195
... Procrustes: rmse 0.04649384 max resid 0.2115397
Run 45 stress 0.1502646
Run 46 stress 0.1518985
Run 47 stress 0.1497692
... Procrustes: rmse 0.03326784 max resid 0.1546212
Run 48 stress 0.1509269
Run 49 stress 0.1502629
Run 50 stress 0.1533299
Run 51 stress 0.1544682
Run 52 stress 0.1541878
Run 53 stress 0.1541584
Run 54 stress 0.1501233
Run 55 stress 0.1501996
Run 56 stress 0.1492999
... New best solution
... Procrustes: rmse 0.001417234 max resid 0.006130985
... Similar to previous best
*** Solution reached
T2_rare<-ggordiplots::gg_ordiplot(ord = ord_t2_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T2")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T2_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_rare_T2.pdf")
Saving 7.29 x 4.51 in image
ord_t3_rare<-ordinate(physeq = subset_samples(ps_rare, Time=="T3"), method = "NMDS", distance = "bray", k=3, trymax= 100)
Square root transformation
Wisconsin double standardization
Run 0 stress 0.1249587
Run 1 stress 0.1256869
Run 2 stress 0.1257678
Run 3 stress 0.127937
Run 4 stress 0.1278503
Run 5 stress 0.124948
... New best solution
... Procrustes: rmse 0.00253675 max resid 0.00957671
... Similar to previous best
Run 6 stress 0.1258153
Run 7 stress 0.1256208
Run 8 stress 0.1257204
Run 9 stress 0.1255412
Run 10 stress 0.1250279
... Procrustes: rmse 0.006195166 max resid 0.03810374
Run 11 stress 0.1250437
... Procrustes: rmse 0.007767753 max resid 0.03872528
Run 12 stress 0.1258564
Run 13 stress 0.1249904
... Procrustes: rmse 0.01078475 max resid 0.03446339
Run 14 stress 0.1256788
Run 15 stress 0.1251965
... Procrustes: rmse 0.01742177 max resid 0.0711806
Run 16 stress 0.1258553
Run 17 stress 0.1249628
... Procrustes: rmse 0.009053109 max resid 0.02873361
Run 18 stress 0.1256671
Run 19 stress 0.1249994
... Procrustes: rmse 0.0125122 max resid 0.04112397
Run 20 stress 0.1257192
*** Solution reached
T3_rare<-ggordiplots::gg_ordiplot(ord = ord_t3_rare, groups = data.frame(sample_data(subset_samples(ps_rare, Time == "T3")))$Herbicide, choices = c(1, 2), kind = c("se"), conf = 0.95, show.groups = "all", ellipse = TRUE, label = FALSE, hull = FALSE, spiders = FALSE, plot = TRUE, pt.size = 1)
T3_rare$plot + theme_classic()
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_ordination_rare_T3.pdf")
Saving 7.29 x 4.51 in image
box and whisker plots of distance within group distances
library(micrUBIfuns)
beta_boxplot(physeq = subset_samples(ps_rare, Time=="T1"), method = "bray", group = "Herbicide")
$data
$plot
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T1_rare_withingroup_beta.pdf")
Saving 7.29 x 4.51 in image
beta_boxplot(physeq = subset_samples(ps_rare, Time=="T2"), method = "bray", group = "Herbicide")
$data
$plot
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T2_rare_withingroup_beta.pdf")
Saving 7.29 x 4.51 in image
beta_boxplot(physeq = subset_samples(ps_rare, Time=="T3"), method = "bray", group = "Herbicide")
$data
$plot
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T3_rare_withingroup_beta.pdf")
Saving 7.29 x 4.51 in image
treatment to control
plotDistances = function(p, m, s, d) {
# calc distances
wu = phyloseq::distance(p, m)
wu.m = melt(as.matrix(wu))
# remove self-comparisons
wu.m = wu.m %>%
filter(as.character(Var1) != as.character(Var2)) %>%
mutate_if(is.factor,as.character)
# get sample data (S4 error OK and expected)
sd = data.frame(sample_data(p)) %>%
select(s, d) %>%
mutate_if(is.factor,as.character)
sd$Herbicide <- factor(sd$Herbicide, levels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
# combined distances with sample data
colnames(sd) = c("Var1", "Type1")
wu.sd = left_join(wu.m, sd, by = "Var1")
colnames(sd) = c("Var2", "Type2")
wu.sd = left_join(wu.sd, sd, by = "Var2")
#remove this line to plot all comparisons.
wu.sd = wu.sd %>% filter(Type1 == "Hand" | Type1 == "Non-Treated")
# plot
ggplot(wu.sd, aes(x = Type2, y = value)) +
theme_bw() +
geom_point() +
geom_boxplot(aes(color = ifelse(Type1 == Type2, "red", "black"))) +
scale_color_identity() +
facet_wrap(~ Type1, scales = "free_x") +
theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) +
ggtitle(paste0("Distance Metric = ", m))
}
a<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T1"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
a <- a + ggtitle("Time 1 Bray-Curtis Dissimlarities")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T1_rare_allgroup_beta.pdf")
Saving 7 x 7 in image
b<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T2"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
b <-b + ggtitle("Time 2 Bray-Curtis Dissimlarities")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T2_rare_allgroup_beta.pdf")
Saving 7 x 7 in image
c<-plotDistances(p = subset_samples(physeq= ps_rare, Time=="T3"), m = "bray", s = "Barcode_ID_G", d = "Herbicide")
c<- c + ggtitle("Time 3 Bray-Curtis Dissimlarities")
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_T3_rare_allgroup_beta.pdf")
Saving 7 x 7 in image
library(ggpubr)
ggarrange(a, b, c, ncol = 1)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/16S_combined_rare_allgroup_beta.pdf", width = 7, height = 10)
Taxon abundance bar plot
#create super long color vector
col_vector <- c("#000000", "#FFFF00", "#1CE6FF", "#FF34FF", "#FF4A46", "#008941", "#006FA6", "#A30059",
"#FFDBE5", "#7A4900", "#0000A6", "#63FFAC", "#B79762", "#004D43", "#8FB0FF", "#997D87",
"#5A0007", "#809693", "#FEFFE6", "#1B4400", "#4FC601", "#3B5DFF", "#4A3B53", "#FF2F80",
"#61615A", "#BA0900", "#6B7900", "#00C2A0", "#FFAA92", "#FF90C9", "#B903AA", "#D16100",
"#DDEFFF", "#000035", "#7B4F4B", "#A1C299", "#300018", "#0AA6D8", "#013349", "#00846F",
"#372101", "#FFB500", "#C2FFED", "#A079BF", "#CC0744", "#C0B9B2", "#C2FF99", "#001E09",
"#00489C", "#6F0062", "#0CBD66", "#EEC3FF", "#456D75", "#B77B68", "#7A87A1", "#788D66",
"#885578", "#FAD09F", "#FF8A9A", "#D157A0", "#BEC459", "#456648", "#0086ED", "#886F4C",
"#34362D", "#B4A8BD", "#00A6AA", "#452C2C", "#636375", "#A3C8C9", "#FF913F", "#938A81",
"#575329", "#00FECF", "#B05B6F", "#8CD0FF", "#3B9700", "#04F757", "#C8A1A1", "#1E6E00",
"#7900D7", "#A77500", "#6367A9", "#A05837", "#6B002C", "#772600", "#D790FF", "#9B9700",
"#549E79", "#FFF69F", "#201625", "#72418F", "#BC23FF", "#99ADC0", "#3A2465", "#922329",
"#5B4534", "#FDE8DC", "#404E55", "#0089A3", "#CB7E98", "#A4E804", "#324E72", "#6A3A4C",
"#83AB58", "#001C1E", "#D1F7CE", "#004B28", "#C8D0F6", "#A3A489", "#806C66", "#222800",
"#BF5650", "#E83000", "#66796D", "#DA007C", "#FF1A59", "#8ADBB4", "#1E0200", "#5B4E51",
"#C895C5", "#320033", "#FF6832", "#66E1D3", "#CFCDAC", "#D0AC94", "#7ED379", "#012C58",
"#7A7BFF", "#D68E01", "#353339", "#78AFA1", "#FEB2C6", "#75797C", "#837393", "#943A4D",
"#B5F4FF", "#D2DCD5", "#9556BD", "#6A714A", "#001325", "#02525F", "#0AA3F7", "#E98176",
"#DBD5DD", "#5EBCD1", "#3D4F44", "#7E6405", "#02684E", "#962B75", "#8D8546", "#9695C5",
"#E773CE", "#D86A78", "#3E89BE", "#CA834E", "#518A87", "#5B113C", "#55813B", "#E704C4",
"#00005F", "#A97399", "#4B8160", "#59738A", "#FF5DA7", "#F7C9BF", "#643127", "#513A01",
"#6B94AA", "#51A058", "#A45B02", "#1D1702", "#E20027", "#E7AB63", "#4C6001", "#9C6966",
"#64547B", "#97979E", "#006A66", "#391406", "#F4D749", "#0045D2", "#006C31", "#DDB6D0",
"#7C6571", "#9FB2A4", "#00D891", "#15A08A", "#BC65E9", "#FFFFFE", "#C6DC99", "#203B3C",
"#671190", "#6B3A64", "#F5E1FF", "#FFA0F2", "#CCAA35", "#374527", "#8BB400", "#797868",
"#C6005A", "#3B000A", "#C86240", "#29607C", "#402334", "#7D5A44", "#CCB87C", "#B88183",
"#AA5199", "#B5D6C3", "#A38469", "#9F94F0", "#A74571", "#B894A6", "#71BB8C", "#00B433",
"#789EC9", "#6D80BA", "#953F00", "#5EFF03", "#E4FFFC", "#1BE177", "#BCB1E5", "#76912F",
"#003109", "#0060CD", "#D20096", "#895563", "#29201D", "#5B3213", "#A76F42", "#89412E",
"#1A3A2A", "#494B5A", "#A88C85", "#F4ABAA", "#A3F3AB", "#00C6C8", "#EA8B66", "#958A9F",
"#BDC9D2", "#9FA064", "#BE4700", "#658188", "#83A485", "#453C23", "#47675D", "#3A3F00",
"#061203", "#DFFB71", "#868E7E", "#98D058", "#6C8F7D", "#D7BFC2", "#3C3E6E", "#D83D66",
"#2F5D9B", "#6C5E46", "#D25B88", "#5B656C", "#00B57F", "#545C46", "#866097", "#365D25",
"#252F99", "#00CCFF", "#674E60", "#FC009C", "#92896B")
phylumGlommed <- tax_glom(ps_rare, "Phylum")
#t1
phylumGlommed_herb_t1 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T1"), group = "Herbicide")
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
phylumGlommed_herb_t1 <- transform_sample_counts(phylumGlommed_herb_t1, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t1)$Herbicide <- factor(sample_data(phylumGlommed_herb_t1)$Herbicide, levels = c(1, 2, 3, 4, 5),
labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
plot_bar(phylumGlommed_herb_t1, x = "Herbicide", fill = "Phylum") + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) +
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_PT1/Figures/Taxon_barplot_t1.pdf")
Saving 7.29 x 4.51 in image
#t2
phylumGlommed_herb_t2 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T2"), group = "Herbicide")
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
phylumGlommed_herb_t2 <- transform_sample_counts(phylumGlommed_herb_t2, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t2)$Herbicide <- factor(sample_data(phylumGlommed_herb_t2)$Herbicide, levels = c(1, 2, 3, 4, 5),
labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
plot_bar(phylumGlommed_herb_t2, x = "Herbicide", fill = "Phylum") + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) +
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_Pt1/Figures/Taxon_barplot_t2.pdf")
Saving 7.29 x 4.51 in image
#t3
phylumGlommed_herb_t3 <- merge_samples(subset_samples(physeq= phylumGlommed, Time=="T3"), group = "Herbicide")
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
Warning in asMethod(object) : NAs introduced by coercion
phylumGlommed_herb_t3 <- transform_sample_counts(phylumGlommed_herb_t3, function(OTU) OTU/sum(OTU))
sample_data(phylumGlommed_herb_t3)$Herbicide <- factor(sample_data(phylumGlommed_herb_t3)$Herbicide, levels = c(1, 2, 3, 4, 5),
labels = c("Non-Treated", "Hand", "Aatrex", "Clarity", "Roundup Powermax"))
plot_bar(phylumGlommed_herb_t3, x = "Herbicide", fill = "Phylum") + theme_classic() + ggtitle("Proportional Taxon Abundances Time 1") +
theme(legend.position="bottom") + guides(fill=guide_legend(nrow=6)) + geom_bar(stat="identity") + theme(axis.text.x=element_text(angle = 45, hjust = 1, size = 5)) +
scale_fill_manual(values = col_vector)
ggsave("/Users/gordoncuster/Desktop/Git_Projects/Herbicide_Microbes_Pt1/Figures/Taxon_barplot_t3.pdf")
Saving 7.29 x 4.51 in image